home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / gle / util / contour / main.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-29  |  10.8 KB  |  454 lines
  1. /*
  2.     Contouring program, creates gle files.
  3. */
  4.  
  5. #include <ctype.h>
  6. #include <stdio.h>
  7. #include <stdlib.h>
  8. #include <math.h>
  9. #include <string.h>
  10. #define false (0)
  11. #define true (!false)
  12.  
  13. char *clab[100];
  14. float cval[100];
  15. char outlab[80];
  16. char outdat[80];
  17. char outgle[80];
  18. char outkey[80];
  19. char outfile[80];
  20. int smooth=false,smoothsub;
  21. FILE *fp,*flab,*fdat;
  22.  
  23. #ifdef unix 
  24. #define VAXC 1
  25. #endif
  26. #ifdef VAXC
  27. char *strdup(char *s);
  28. #else
  29. #include <alloc.h>
  30. #endif
  31. int gen_next(char *pat, double *v);
  32. int gen_cval(char *cstr, float cval[], int *ncv) ;
  33. FILE *myfopen(char *fname, char *mode);
  34. int read_data(char *fname, float *zmin,float *zmax,int *nx,int *ny);
  35. int docontour(float zz[], long nrz, long nx, long ny
  36.     , float cval[], char *clab[], long ncv, float zmax);
  37.  
  38. int addvect(int m, double x, double y);
  39. FILE *myfopen(char *fname, char *mode)
  40. {
  41.     FILE *f;
  42.     f = fopen(fname,mode);
  43.     if (f==NULL) {
  44.         printf("Unable to open {%s} \n",fname); perror(""); exit(1);
  45.     }
  46.     return f;
  47. }
  48. int dflt(char *s, char *d);
  49. double dfltv( char *d);
  50. double dfltv( char *d)
  51. {
  52.     char buff[80];
  53.     dflt(buff,d);
  54.     return atof(buff);
  55. }
  56. dflt(char *s, char *d)
  57. {
  58.     printf(" [%s] ? ",d); gets(s);
  59.     if (strcmp(s,"")==0) strcpy(s,d);
  60. }
  61. #ifdef VAXC
  62. #define farmalloc malloc
  63. #endif
  64. static double dxmin,dymin,dxmax,dymax;
  65. static float *zz;
  66. static     double nnx,nny;
  67. main()
  68. {
  69.     float zmin = 1e10, zmax = -1e10;
  70.     char infile[80],cdflt[80],cstr[80],lablet[80];
  71.     int nx,ny,i,ncv,dokey;
  72.     int zdim=51;
  73.     char buff[90];
  74.     
  75.     printf("Data file containing Z values");  dflt(infile,"");
  76.     read_data(infile,&zmin,&zmax,&nx,&ny);
  77.     zdim = nx;
  78.     if (dxmin==dxmax) {dxmax = nx; dxmin = 1;}
  79.     if (dymin==dymax) {dymax = ny; dymin = 1;}
  80.  
  81.     sprintf(cdflt,"%g:%g:%g",zmin,zmax,(zmax-zmin)/10.0);
  82.     printf("Specify values to contour at, e.g. 1 2 3 4 or 1:4:.5\n");
  83.     dflt(cstr,cdflt);
  84.     printf("Label contours with letters A...Z");  dflt(lablet,"YES");
  85.     if (toupper(lablet[0]) == 'Y') dokey = true; else dokey = false;
  86.     gen_cval(cstr, cval, &ncv);
  87.     for (i=0; i<ncv; i++) {
  88.       if (dokey) {
  89.         sprintf(buff,"%c",i+'A');
  90.         clab[i] = strdup(buff);
  91.         printf("{%s} ",clab[i]);        
  92.       } else {
  93.         sprintf(buff,"%g",cval[i]);
  94.         clab[i] = strdup(buff);
  95.       }
  96.     }     
  97.     if (toupper(lablet[0]) == 'Y') printf("\n");
  98.     printf("Contour values ");
  99.     for (i=0; i<ncv; i++)    printf("%g ",cval[i]);
  100.     printf("\n");
  101.     printf("Warning: Smoothing is SLOW\n");
  102.     printf("Smoothing factor (points per grid step 3-10, 0=no smoothing) "); dflt(lablet,"0");
  103.     smoothsub = atoi(lablet);
  104.     if (smoothsub>0) smooth = true;
  105.  
  106.     printf("Output file names (E.G. xxx.gle, xxx.dat, xxx.lab, xxx.key");
  107.     dflt(outfile,"xxx");
  108.     
  109.     strcpy(outgle,outfile); strcat(outgle,".gle");
  110.     strcpy(outdat,outfile); strcat(outdat,".dat");
  111.     strcpy(outlab,outfile); strcat(outlab,".lab");
  112.     strcpy(outkey,outfile); strcat(outkey,".key");
  113.  
  114.     fp = myfopen(outgle,"w");
  115.     fprintf(fp,"size 24 18\n");
  116.     fprintf(fp,"begin graph\n");
  117.     fprintf(fp,"    size 24 18\n");
  118.     fprintf(fp,"    xaxis min %g max %g  \n",dxmin,dxmax);
  119.     fprintf(fp,"    yaxis min %g max %g  \n",dymin,dymax);
  120.     fprintf(fp,"    d1 lstyle 1 \n",nx,ny);
  121.     fprintf(fp,"    d1 bigfile %s \n",outdat);
  122.     fprintf(fp,"end graph\n");
  123.     fprintf(fp,"save theorigin\n");
  124.     fprintf(fp,"set hei .3 just cc\n");
  125.     fprintf(fp,"bigfile %s\n",outlab);
  126.     fprintf(fp,"move theorigin\n");
  127.     fprintf(fp,"set hei .3 just left \n");
  128.     if (dokey) fprintf(fp,"bigfile %s\n",outkey);
  129.     fclose(fp);
  130.  
  131.     flab = myfopen(outlab,"w");
  132.     fdat = myfopen(outdat,"w");
  133.     fprintf(flab,"set lwidth .0000001\n");
  134.     nnx = nx; nny = ny;
  135.  
  136.     docontour(zz,zdim,nx,ny,cval,clab,ncv,zmax);
  137.  
  138.     fclose(flab);
  139.     fclose(fdat);
  140.  
  141.     fp = myfopen(outkey,"w");
  142.     fprintf(fp,"rmove .2 .2 \n");
  143.     fprintf(fp,"set lwidth 0\n");
  144.     fprintf(fp,"save keystart \n");
  145.     fprintf(fp,"begin box fill white add .2 \n");
  146.     for (i=ncv-1; i>=0; i--) {
  147.         fprintf(fp,"text %s\n",clab[i]);
  148.         fprintf(fp,"rmove .6 0 \n");
  149.         fprintf(fp,"text %g\n",cval[i]);
  150.         fprintf(fp,"rmove -.6 .5\n");
  151.     }
  152.     fprintf(fp,"end box \n");
  153.     fprintf(fp,"move keystart \n");
  154.     for (i=ncv-1; i>=0; i--) {
  155.         fprintf(fp,"text %s\n",clab[i]);
  156.         fprintf(fp,"rmove .6 0 \n");
  157.         fprintf(fp,"text %g\n",cval[i]);
  158.         fprintf(fp,"rmove -.6 .5\n");
  159.     }
  160.     fprintf(fp,"move keystart \n");
  161.     fclose(fp);
  162. }
  163. docontour(float zz[], long nrz, long nx, long ny
  164.     , float cval[], char *clab[], long ncv, float zmax)
  165. {
  166.     extern int draw_(float *x, float *y, long *iflag), plot_(float *, float *, int *);
  167.     static long *work, i, j;
  168.     extern  gcontr_(float *, long *, long *, long *, float *, long *, float *, long *, 
  169.        int draw_(float *x, float *y, long *iflag) );
  170. /*     dimension of work is large enough to contain     */
  171. /*     2*(dimension of c)*(total dimension of z) useful bits.  */
  172.     work = farmalloc(( (long) sizeof(long) * 2 * ncv * nx * ny)/31 + 10);
  173.     if (work==NULL) {
  174.         printf("Unable to allocate storage for work array\n");
  175.         exit(1);
  176.     }
  177.     zmax = zmax + 100; /* no clipping */
  178.         gcontr_(zz, &nrz, &nx, &ny, cval, &ncv, &zmax, work, draw_);
  179. }
  180.  
  181. /*#define sx(v) (dxmin + (dxmax-dxmin)*(v)/nnx) */
  182. #define sx(v) (dxmin + (dxmax-dxmin)*(v-1)/(nnx-1))
  183. #define sy(v) (dymin + (dymax-dymin)*(v-1)/(nny-1))
  184. int draw_(float *x, float *y, long *iflag)
  185. {
  186.     static float xcur,ycur;
  187.     static int ih, il;
  188.  
  189.     ih = *iflag / 10;
  190.     il = *iflag - ih * 10;
  191.     switch (il) {
  192.       case 6:   /* Get current point */
  193.         *x = xcur;
  194.             *y = ycur;
  195.             break;
  196.       case 2: /* start contour at bondary */
  197.       case 3: /* start contour not at bondary */
  198.         if (smooth) addvect(1,sx(*x),sy(*y));
  199.         else {
  200.             fprintf(fdat,"* * \n",*x,*y);
  201.             fprintf(fdat,"%g %g \n",sx(*x),sy(*y));
  202.         }
  203.         fprintf(flab,"amove xg(%g) yg(%g) \n",sx(*x),sy(*y));
  204.         fprintf(flab,"begin box fill white add .07 \n");
  205.         fprintf(flab," text %s \n",clab[ih-1]);
  206.         fprintf(flab,"end box \n");
  207.         fprintf(flab,"text %s \n",clab[ih-1]);
  208.         break;
  209.       case 1: /* continue a contour */
  210.         if (smooth) addvect(2,sx(*x),sy(*y));
  211.         else fprintf(fdat,"%g %g \n",sx(*x),sy(*y));
  212.         break;
  213.       case 4: /* finish contour at a boundary */
  214.         if (smooth) addvect(4,sx(*x),sy(*y));
  215.         else fprintf(fdat,"%g %g \n",sx(*x),sy(*y));
  216.         break;
  217.       case 5: /* finish close contour */
  218.         if (smooth) addvect(3,sx(*x),sy(*y));
  219.         else fprintf(fdat,"%g %g \n",sx(*x),sy(*y));
  220.         break;
  221.     }
  222.     xcur = *x;
  223.     ycur = *y;
  224. } 
  225. gen_cval(char *cstr, float cval[], int *ncv)
  226. {
  227.     double v;
  228.     int j;
  229.     *ncv = 0;
  230.     gen_next(cstr,&v);
  231.     for (;gen_next(NULL,&v);) {
  232.         cval[*ncv] = v;
  233.         (*ncv)++;
  234.     }
  235. }
  236. gen_next(char *pat, double *v)
  237. {
  238.     static char *s,p[200];
  239.     static double v1,v2,v3;
  240.     static int instep;
  241.     char *c1,*c2;
  242.     if (pat!=NULL) {
  243.         strcpy(p,pat);
  244.         s = strtok(p,", ");
  245.         return true;
  246.     }
  247.     if (instep) {
  248.            v1 += v3;
  249.            *v = v1;
  250.            if (v1<=(v2+v3/10)) return true;
  251.            else instep = false;
  252.     }
  253.     if (s==NULL) return false;
  254.     c2 = c1 = strchr(s,':');
  255.     if (c1!=NULL) c2 = strchr(c1+1,':');
  256.  
  257.     if (c1==NULL) {
  258.         *v = atof(s);
  259.         s = strtok(NULL,", ");
  260.         return true;
  261.     }
  262.     v1 = atof(s);
  263.     v2 = atof(c1+1);
  264.     if (c2==NULL) v3 = 1; else v3 = atof(c2+1);
  265.     *v = v1;
  266.     instep = true;
  267.     s = strtok(NULL,", ");
  268.     if (v1>v2) return false;
  269.     return true;
  270. }
  271. static char buff[2032];
  272. FILE *myfopen(char *fname, char *mode);
  273. int alloc_zdata(int nx, int ny);
  274. alloc_zdata(int nx, int ny)
  275. {
  276.     if (zz!=NULL) free(zz);
  277.     zz = farmalloc((long) nx * ((long) ny+1) * sizeof(float));
  278.     if (zz==NULL) {
  279.         printf("Unable to allocate enough memory for datafile\n");
  280.         return true;
  281.     }
  282.     return false;
  283. }
  284. double getkeyval(char *buff,char *k);
  285. double getkeyval(char *buff,char *k)
  286. {
  287.     char *s;
  288.     s = strstr(buff,k);
  289.     if (s!=NULL) return atof(s+strlen(k));
  290.     return 0.0;
  291. }
  292. FILE *df;
  293. read_data(char *fname, float *zmin,float *zmax,int *nx,int *ny)
  294. {    
  295.     double v;
  296.     long x,y;
  297.     int c,b;
  298.     char *s;
  299.     int i;
  300.     x = y = 0;
  301.     *nx = 0; *ny = 0;
  302.  
  303.     df = myfopen(fname,"r");
  304.     if (df==NULL) {*nx = 0; *ny = 0; return;}
  305.     for (;!feof(df);) {
  306.       if (fgets(buff,2000,df)!=NULL) {
  307.         if (*nx==0) {
  308.             strupr(buff);
  309.             *nx  = getkeyval(buff,"NX");
  310.             *ny  = getkeyval(buff,"NY");
  311.             dxmin = getkeyval(buff,"XMIN");
  312.             dymin = getkeyval(buff,"YMIN");
  313.             dxmax = getkeyval(buff,"XMAX");
  314.             dymax = getkeyval(buff,"YMAX");
  315.             if (*nx==0 || *ny==0) {
  316.                 printf("Expecting to find  ! nx 11 ny 11 xmin 0 xmax 10 ymin 0 ymax 10 \n");
  317.                 printf("on first line of data file,  Please enter this data \n");
  318.                 printf("Number of columns (nx) ");  *nx = dfltv("");
  319.                 printf("Number of rows (ny) ");  *ny = dfltv("100");
  320.                 printf("Min x ");  dxmin = dfltv("0");
  321.                 printf("Max x ");  dxmax = dfltv("1");
  322.                 printf("Min y ");  dymin = dfltv("0");
  323.                 printf("Max y ");  dymax = dfltv("1");
  324.             } else {
  325.                 fgets(buff,2000,df);
  326.             }
  327.             if (alloc_zdata(*nx,*ny)) return;
  328.         }
  329. check_again:
  330.         b = strlen(buff);
  331.         c = buff[b-1];
  332.         if (strchr(" \n\t",c)==NULL) { /* we're halfway thru a number */
  333.             buff[b] = getc(df);
  334.             buff[b+1] = 0;
  335.             goto check_again;
  336.         }
  337.         s = strchr(buff,'!');
  338.         if (s!=NULL) *s = 0;
  339.         s = strtok(buff," \t\n,");
  340.         for (;s!=NULL;) {
  341.             v = atof(s);
  342.             if (isdigit(*s) || *s=='-' || *s=='+' || *s=='.') {
  343.                 if (x>= *nx) {
  344.                     x = 0; y++;
  345.                 }
  346.                 if (y>= *ny) {
  347.                     printf("Too much data in data file %ld %d \n",y,*ny);
  348.                     return;
  349.                 }
  350.                 if (v < *zmin) *zmin = v;
  351.                 if (v > *zmax) *zmax = v;
  352.                 zz[x + y * *nx] = v;        
  353.                 x++;
  354.             } else printf("Not a number {%s} \n",s);
  355.             s = strtok(NULL," \t\n,");
  356.         }    
  357.       }
  358.     }
  359.     fclose(df);
  360.     *ny = y+1;
  361.     return;
  362. abort_data:
  363.     fclose(df);
  364.     *nx = 0; *ny = 0;
  365.     return;
  366.  
  367. }
  368.  
  369. #ifndef __TURBOC__
  370. char *strdup(char *s)
  371. {
  372.     char *ss;
  373.     ss = malloc(strlen(s)+1);
  374.     strcpy(ss,s);
  375.     return ss;
  376. }
  377. strupr(char *s)
  378. {
  379.     for (;*s!=0;s++) *s = toupper(*s);
  380. }
  381. #endif
  382. plot_(float *a, float *b, int *c){}
  383.  
  384.  
  385. int nsm,nout_alloc;
  386. float xin[400],yin[400],*xout,*yout;
  387. addvect(int m, double x, double y)
  388. {
  389.     long mode,ninp,nsub,nout;
  390.     int i;
  391.     int snipends=false;
  392.     if (m==1) nsm = 0;
  393.     if (nsm>0) if (xin[nsm-1]==x && yin[nsm-1]==y) {
  394.         return;
  395.     } 
  396.     xin[nsm] = x;
  397.     yin[nsm++] = y;
  398.     if (nsm>299) {
  399.         printf("Too many points to smooth, don't use smooth option\n");
  400.         exit(1);
  401.     }
  402.     /* Connect the ends by adding some points if contour is right sort */
  403.     /* 4=boundary, 3=closed shape */
  404.     if (m==3 || m==4) { 
  405.         if (nsm<2) {
  406.             for (i=0;i<nsm;i++) {
  407.                 fprintf(fdat,"%g %g \n",xin[i],yin[i]);
  408.             }
  409.             nsm = 0;
  410.             return;
  411.         }
  412.         if (m==3) {
  413.             snipends = true;
  414.             for (i=nsm; i>0; i--) { 
  415.                 xin[i] = xin[i-1];
  416.                 yin[i] = yin[i-1]; 
  417.             }
  418.             xin[0] = xin[nsm-1];
  419.             yin[0] = yin[nsm-1];
  420.             xin[nsm+1] = xin[2];
  421.             yin[nsm+1] = yin[2];
  422.             nsm += 2;
  423.             /*printf("\n");
  424.             for (i=0; i<nsm; i++) { 
  425.                 printf("%g %g \n",xin[i],yin[i]);
  426.             }*/
  427.         }
  428.  
  429.         ninp = nsm;        /* Number of input points */
  430.         mode = 2;         /* multi valued */
  431.         nsub = smoothsub;    /* user parameter */
  432.         if (nsub<3) nsub = 3;
  433.         nout = (ninp-1)*nsub+1;
  434.         if (nout>nout_alloc) {
  435.             xout = malloc(nout*sizeof(*xout));
  436.             yout = malloc(nout*sizeof(*xout));
  437.             nout_alloc = nout;
  438.         }
  439.         glefitcf_(&mode,xin,yin,&ninp,&nsub,xout,yout,&nout);
  440.         nsm = 0;
  441.  
  442.         fprintf(fdat,"* * \n");
  443.         if (snipends) {
  444.          for (i=nsub;i<nout-nsub;i++) {
  445.             fprintf(fdat,"%g %g \n",xout[i],yout[i]);
  446.          }
  447.         } else {
  448.          for (i=0;i<nout;i++) {
  449.             fprintf(fdat,"%g %g \n",xout[i],yout[i]);
  450.          }
  451.         }
  452.     }
  453. }
  454.